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Abstract 

The shock-layer flowfield is obtained with or 
without viscous and heat-conducting dissipations 
from the conservative laws of fluid dynamics 
equations using a shock-fitting implicit finite- 
difference technique. The governing equations are 
cast in curvilinear-orthogonal coordinates and 
transformed to the domain between the shock and 
the body. Another set of equations is used for 
the singular coordinate axis, which, together with 
a cone generator away from the stagnation point, 
encloses the computation domain. After initializ- 
ing the flow variables on a prescribed grid, a 
time-dependent alternating direction implicit fac- 
torization technique is applied to integrate the 
equations with local-time increment until a steady 
solution is reached. The shock location is up- 
dated after the flowfield computation, but the 
wall conditions are implemented into the implicit 
procedure. Since primitive variables and few 
metrics are used, the numerical formulation is 
simple and the core requirement is not stringent. 
Innovative procedures have been introduced to 
define the initial flowfield, to treat both 
perfect and equilibrium gases, to advance the 
solution on a coarse- to-fine grid sequence, and to 
start viscous flow computations from their corre- 
sponding inviscid solutions. The code has proven 
capabilities for a wide range of free-stream 
conditions and body configurations. Among the 
examples shown are the Space Shuttle Orbiter 
equilibrium flow case at Mach 22 and an angle of 
attack of 40.8®, and an aerobraking orbital trans- 
fer vehicle perfect-gas case having a 60® cone and 
sonic shoulder at Mach 34.8. These results are 
obtained from a grid no greater than 28 by 18 by 7 
and converged within 300 integration steps. They 
are of sufficient accuracy to start parabolized 
Navier-Stokes or Euler calculations beyond the 
nose region, to compare with flight and wind- 
tunnel data, and to evaluate conceptual designs of 
reentry spacecraft. 

1 . Introduction 

Although the solutions of blunt-body flow 
have been pursued ever since the beginning of 
supersonic flight, they received most attention 
later in the design of spacecraft to alleviate the 
aerodynamic heating during entry. If the viscous 
effect is confined to a thin boundary layer adja- 
cent to the body, the inviscid assumption can be 
invoked to simplify the complexity of the problem. 
Even so, the governing equations are nonlinear and 
of mixed characteristics bounded by the bow shock; 
hence, only numerical methods are applicable. The 
classical approaches for inviscid flow computation 
developed before I960 have been summarized in the 
book by Hayes and Probstein,0) wherein three- 
dimensional (30) problems were barely touched 
because of the limited computer resources avail- 
able then. The new era of blunt-body flow com- 
putation began with the pioneering work of 
Moretti,(2) who introduced a shock-fitting proce- 
dure to the well-formulated time-dependent finite- 
difference approach. The problem has become a 


time-dependent, or iterative, solution starting 
from a reasonable guess of shock location and 
shock-layer variables, whereby a firm theory of 
initial-value problems can be used. Moretti and 
Bleich's(2) three-dimensional code was made avail- 
able in 1967. Since then, it was improved and 
adopted by many organizations for their specific 
applications. The most important ones, probably, 
are the finite-volume code by Rizzi and Inouye,(o) 
the A-scheme code by Hall,(^) and the split- 
coefficient code by Daywitt.(5) Whereas smooth 
bodies are the objective aimed at in Refs. 2 and 
3, indented nosetips are of primary interest as 
investigated in Refs. 4 and 5. These Euler codes 
are based on explicit technique and confined to 
the shock layer, a common feature shared by many 
other versions developed in the 1970 decade. 

The development of three-dimensional Navier- 
Stokes (30 NS) codes has progressed relatively 
slowly but still shown impressive accomplishment. 
The first attempt was made by the author (6) in 
1974, when he used MacCormack*s two-step scheme 
and a finite-difference procedure to calculate the 
shock and body variables as opposed to using 
Moretti *s characteristic procedure. Since dense 
grid points used to resolve the diffusive fluxes 
near the wall penalized the rate of convergence, 
results were unsatisfactory even at the cost of 
several hundred hours of IBM 360-67 computer 
time. (7) Obviously, after demonstrating the 
feasibility of NS methodology, more powerful 
numerical techniques and computers are needed to 
perform any computation for practical purposes. 

As soon as Beam and Warming's alternating direc- 
tion implicit (ADI) factorization technique(Q) 
became fully developed, it was adapted and 
modified to be suitable for the peculiarity of 
blunt-body problems. The computation time was 
reduced to a mere 10 hr of Univac 1182 computer 
because of the use of less restricted time-step 
size. The implicit code was used successfully to 
solve some very difficult problems for the Space 
Shuttle Orbiter at high angles of attack, and 
preliminary results have been reported in 1981.(9) 
During this period, Kutler et al.OO) developed a 
3D shock-fitting NS code on the basis of Pulliam 
and Steger'sHI) generalized-coordinate NS code. 
Since the motivation of Kutler 's work came from 
the need to predict aerodynamic characteristics 
for an indented nose tip at relatively low angles 
of attack, the code is known to have had stability 
difficulties and failed to converge for flow 
incidence angles of greater than 30®. 

The purpose of this paper is to present the 
detailed formulation, the method of computation, 
and the verification sample cases of the 3D NS 
code, which has undergone continued refinement 
since 1981. A substantial acceleration of conver- 
gence rate has been realized by means of implicit 
damping and of local-time step at a constant 
Courant number (CN). The accuracy of equilibrium- 
air solution is enhanced by using Tannehill’s 
routines for the eauation of state and the trans- 
port coefficients.v 12, 13) Additional schemes to 
sequence the computations from coarse to fine 


1 



grids are also found useful to shorten the number 
of iteration cycles. The performance of the 3D NS 
code is much superior to that of its original 
version, the capabilities of considering different 
body geometries have been broadened, and strong 
flow asymmetry can be easily handled. 

The code is intended to meet two main re- 
quirements: prediction of wall heat-flux distri- 

bution and definition of flow variables on an 
axial plane in the supersonic region. The NS 
methodology is applicable to blunt reentry space- 
craft at high altitude where the low-density and 
low-Reynolds-number conditions invalidate the thin 
boundary-layer (BL) assumption for a decoupled 
Euler and BL approach. The NS code is also needed 
for determining the flowfield around a highly 
blunt body truncated by a sonic shoulder as the 
pressure and boundary layer are strongly influ- 
enced by the abrupt change of body contour. Gen- 
eration of flow variables to start the downstream 
computation by a space-marching technique can be 
conventionally done with an axisymmetric two- 
dimensional (2D) code for a spherical nose. How- 
ever, for nonspherical nose shapes and high angles 
of attack, the 3D NS code should be used. 

The formulation and methods of calculation 
are tailored to satisfy the goals and selected for 
use on scalar computers (Univac 1182) having 3 
million instructions per second (MIPS) and 262 000 
words of main core. In view of the body config- 
urations illustrated in Fig. 1 and the absence of 
embedded shocks, the conventional, weakly conser- 
vative NS equations are cast in curvilinear- 
orthogonal coordinates and converted to the 
spherical-polar frame by two metrics. The singu- 
lar axis is governed by a different set of NS 
equations in order to eliminate solution errors 
introduced by other means of solution. 
straightforward central-difference quotient is 
used to approximate the derivatives of the con- 
vective and diffusive flux terms. The numerical 
solution is time-dependent starting from an 
approximate shock and initial assumption of flow 
variables. The shock is fitted explicitly, but 
the wall and the shock layer itself are updated by 
the ADI factorization technique, which, in turn, 
leads to three systems of linear algebraic equa- 
tions of block-matrix structure. Primitive 
variables are sought for by the implicit solution 
to lessen the amount of arithmetic in evaluating 
Jacobians. The Euler solution is often obtained 
to initialize the NS computation, although, in 
some cases, the direct application of coarse-to- 
fine grid sequence is more cost-effective for NS 
computations. Either is needed to reduce the 
computation cost. 

The remaining portion of the paper is divided 
into four sections and an appendix. The vector 
notation of coordinate systems is defined immedi- 
ately following this section. Duplication of gov- 
erning equations already presented in earlier 
publications (6, 9) cannot be totally avoided for 
the completeness of this final report. In the 
third section, the well-known factorization tech- 
nique is briefly outlined with emphasis placed on 
the Jacobian matrices and damping terms which are 
different from other work. (8> 10, 1 1 ) xhe fourth 
section is intended to summarize the capabilities 


and limitations of the 3D NS code and includes 
remarks on proper use of the code. The discussion 
of results in sec. 5 begins with the explanation 
of important parameters and their effects on the 
solution, then ends with sample results demon- 
strating the code usage. Appendix A contains 
detailed derivation of the governing equations on 
the singular coordinate axis. 

2. Theoretical Formulation 

2.1 Geometric Consideration and Coordinate 
Systems 

Because of the commonness of blunt entry 
bodies, a flowfield solution is often needed to 
provide the heat-flux distributions over a non- 
spherical nose, wide-angle body with sonic 
shoulder and to determine the flow variables on an 
axial plane for computations downstream. Figure 1 
depicts the three basic configurations differing 
in the orientation of the coordinate system and 
the wind direction. The flow incidence angle a 
and the reference angle have either a positive 
value or a zero as shown, or both can have posi- 
tive values for additional convenience to define 
the region of interest. For most cases, there is 
a plane of symmetry so that the flowfield is only 
needed to the right side of the wind vector. 

The computation region extended from the body to 
the shock and enclosed by the Z-axis, and the 
cone generator with various vertex angles are 
seen in Fig. 1. As has been done in practice, a 
spherical-polar system is used to define the 
physical space over the frontal portion of the 
blunt body, then a cylindrical-polar system is 
used downstream of the interface shared by both 
coordinate systems. ( 6 , 9 ) If a more complicated 
body such as an indented nose is to be considered, 
conformal mappings of meridional plane may be used 
to generate a family of coordinate lines less 
oblique to the body wall than those in the 
spherical system. Fortunately, the lack of con- 
trol of the location and the orientation of the 
coordinate lines at the boundary does not pose a 
serious limitation to the computation of smooth 
bodies. Furthermore, the present formulation 
requires only two metrics in contrast to nine 
metrics and a Jacobian associated with the 
generalized-coordinate formulation. Consequently, 
the arithmetic is comparatively reduced. 

Three additional coordinate systems are used 
to complement the physical spherical system for 
the purpose of interrelating the velocity vector 
components. They are referred to as wind, body, 
and intrinsic frames and shown in Fig. 2. The 
unit vectors are defined in Table 1 . 

Let curvilinear coordinates (^,n>0 be ^ = n - 
0, q = r, and < = 4>. Then, the unit vectors 
(iyjfk) and UyJyK) shown in Fig. 2 are related by 


df ^ I dX J dY '^KdZ^id^-¥jdr\’¥kdl, 


and 


X = q sin ^ cos y = q sin ^ sin Z = -q cos ^ 
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The relationships are given in Table 1. If the 
body surface is represented by/ = n-B{^,0 =0, 
then its outward normal gradient is V/ = y - iB^/B 
- kB(^/(B sin 0. After normalization of V/*, n = 
V//|v/|, the tangent and binomial vectors, also 
shown in Fig. 1, are obtained from t = n X k and 6 
^ t X n. The intrinsic frame on the shock ^ = n - 
= 0 is defined likewise, except that n is 
inward. Note that all unit vectors follow the 
right-hand rule. 

The body configuration is prescribed analyt- 
ically in terms of a series of conics or their 
degenerated form assuming symmetry with respect to 
the Z-axis, or on a number of X-Y planes away 
from a spherical cap. The cylindrical (r,4>,Z) 
coordinates are used to define the body, then the 
independent variables are changed to 

2.2 Equations in Computation Space 

The conservative laws of a single-component, 
compressible, viscous flow are given in vector 
notation as follows. 


where 

u = u , f = (2 - u*/u) f + A^-uu*, g = y^h^g 
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( 1 ) 


u* = u(0,ii,0), w = 0 


The derivation of Eq. (2) is given in Appendix A. 


where y ^ ^B) f = A3/, 

g - + yv[h\hzg h = AiA, r = r, = q, 

/13 = n sin and metrics y^y and are held 
nondifferentiable with respect to y. Note that yt 
is absent, even though S = S(^,^,0. The convec- 
tive and diffusive fluxes are of the following 
expression . 


In Eqs, (1) and (2), complete stress tensor 
and heat-flux components are shown, but the order- 
of-magnitude analysis applied to the same equa- 
tions on body intrinsic coordinates indicates that 
and are the only important ones in /, g, h, 
and r. It should be pointed out that in r can- 
not be eliminated at all. 
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In the preceding expressions, p, and e 

are density, velocity components in U,q,0, and 
total energy, respectively. Variables of Eqs. (1) 
and (2) are made dimensionless from 


P = pVp'„ . U = u'/Vp'Jp’l , e = , P = pVp'„ , 


T = r/r„ , l, = . t = t'y/p'Jp'JR^ 

The primed variables indicate dimensional 
quantities, whereas the subscript <» refers to 
ft^ee-stream condition; Rn is the nose radius, and 
PfT are, respectively, pressure and temperature. 

Two components of the stress tensor, ny and 
qiy are defined in orthogonal coordinates as 
follows; 


Equation ( 1 ) is used everywhere except on the 
Z-axis (Fig. 1), wherein w; = 0 and ^ = 0 by defi- 
nition. The exact form of the governing equation 


^ ("« ■ 3 

_ /j_ to 

(y - i)Pj^Re /ij a^ 


u, + /■ + g + r = 0 

t ' X 


( 2 ) 


(3) 


The Stokes assumption has been used to derive 
formulas for the stress tensor. Parameters in Eq. 
(3) are the ratio of specific heats y» Mach number 
M, Reynolds number Re, and Prandtl number Pr. 
Definitions of Re and are given by Re = 
p»Va>Pjv/p and Pr - Cpp/iC; denotes free-stream 
speed, p and K the coefficients of viscosity and 
thermal conductivity, and Cp the specific heat at 
constant pressure. These are dimensional quanti- 
ties. Finally, the thermal equation of state 


P = pRT (4) 


may be solved by means of Newton *s technique. 
p"+‘ = p" + (p^)V-y") 

where 



is required to complete the system of governing 
equations. In this equation, R refers to the gas 
constant. The number of dependent variables is 
limited to six since e = e + 0.5(u2 + u2 + w^) and 
the specific internal energy e can be related to T 
by the caloric equation of state e = CyT, where Gy 
is the specific heat at constant volume, or by e = 
e(p,p) for equilibrium air. 

2.3 Mapping Functions 

Equations (1) and (2) are dedicated to the 
flowfield calculation between the shock and the 
body. This calculation is accomplished by se- 
lecting a coordinate transformation equation such 
as^' = q = (q - B)/{S - P), which implies that 
for a given ^ and q corresponds to equally 
spaced q varying from 0 to 1 . The convenience of 
addressing a radial location within the shock 
layer can be extended to clustered locations of q 
by using any one of the three exponential func- 
tions 


y = (exp (p q ) - l)/(exp (P) - i) 


y = sinh (P q )/sinh p 


^ ln[(P+ q)/(P-q)] 

^ ln[(P + l)/(P-l)l 

Hence, for a given set of y points, uneven distri- 
butions of q points can be defined. However, more 
cluster points toward the body often results in 
fewer points away from the shock. By comparing 
the characteristics due to different exponential 
functions, it is found that the last function from 
TannehillC 14) offers the most desirable clustering 
near y = 1 while maintaining an approximately 
linear relation y « q for y < 0.8. 

As shown in Fig. 3, the degree of clustering 
is determined by the slope of at q = 1, which, 
in turn, is controlled by the parameter p. To 
find p for specified values of y' and q', a non- 
linear equation 


(5a) 

(5b) 

(5c) 


f(P) = y In 


P + i 

P-i 


+ In 


P+ n' 

P- n' 


= 0 


2.4 Boundary Conditions 

Flow variables along q = S(^,0 as given by 
Eqs. (1) and (2) are made to satisfy the Rankine- 
Hugoniot (RH) relations by introducing a local 
shock speed jSt a^ (^>0. Thus, the incoming free- 
stream velocity V* and t^e flow vector immediately 
downstream of the shock V are modified by the 
shock speed in the stationary intrinsic frame. 

Let prime refer to the intrinsic frame. Then 




where = IV^sin (a — -f /fV^cosCa - a^,) 

V = lU + KW = iu + ju + kw 
00 00 00 00 00 


Table 1 can be used to relate ua>t Uoo, and Woo to (/«, 
and Woo. 


An iteration procedure (secant formula) is 
used to find by matching v' with niu + ti 2 V + 
n^w, whereas v' is governed by the RH equations 
shown in the following: 


n(Vl-pV') = 0 

(n-v'y + 1 = p(n-vf + p 

6 (V„ - V') = 0 

HV„-n = 0 

+ 0.5V‘J = h + 0.5V* 

For a given S^, another round of iterations is 
needed to solve Eq. (6). The iterated variable is 
the density ratio across the shock for an ideal- 
gas model; the process converges within three to 
four iterations. For an equilibrium-air model, 
the iterated variable is the enthalpy, which is 
available as a function of A = A(p,p) in a curve- 
fit subroutine in Ref. 12. 
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The wall boundary x\ = imposes differ- 

ent types of conditions depending on whether the 
flow is inviscid or viscous. In accordance with 
the mathematical characteristics of a parabolic 
equation, either pressure or density needs to be 
calculated, whereas other variables may be speci- 
fied. Thus, in the viscous calculation, no-slip 
velocity and isothermal or adiabatic temperature 
are specified on the wall. These conditions are 
simply V = 0, e = e„; or n-Ve a 0. For the hyper- 
bolic equation, only one condition is allowed; 
therefore, n'V = 0 is generally required. In an 
attempt to ensure that the flow variables still 
satisfy conservative laws, th^ velocity components 
parallel to the wall, V'* = and V't = 
are multiplied by a factor V-V/(V'i,2 + then 

converted back to the {^,ri,0 frame by 




The pressure and density as computed from Eqs. (1) 
and (2) remain unchanged. 

The outflow boundary, specified as a cone 
generator with ^ is to be located well 

within the supersonic flow and on the portion of 
the body where a rapid change of flow direction is 
not expected. Inside the viscous layer, when the 
outflow is subsonic, one-point extrapolation is 
used. Otherwise, the flow variables are obtained 
by a linear two-point extrapolation of interior 
results. On the pitch plane K = 0, the flow 
variables are to satisfy k'Vu = 0, where u = py 
p, e, Uy and Vy and to satisfy k V(hVw) = 0. 

The remaining boundary ^ = 0, a singular line 
in physical space, necessitates rigorous analysis. 
The best approach to handling this boundary is to 
solve the governing equations in special form (Eq. 
(2)) on the same basis as for the interior region. 
But Eq. (2) is only to be used once in conjunction 
with flow variables on the pitch plane. Then, 
flow variables on other meridional planes are ob- 
tained from 


p(0,q,0 = p(0,q,0) 
u(0,n»0 = a(0,q,0)cos^ 
u(o,n,0 ~ y(o,ri,o) 

u;(0,ri,0 = -ii(0,ri,0)sm< 
e(0,q,0 = e(0,q,0) 

2.5 Initial Flow Approximation 

The present approach to solving blunt-body 
problems is a special class of time-dependent 
method that retains the unsteady derivative of the 
governing equations but merely simulates pseudo- 
physical phenomena. In a truly numerical simu- 
lation of flow over the body, it is often neces- 
sary to start the time-dependent calculation at 
the instant the body is immersed into the free 
stream, and the calculation continues until the 


steady state is reached. At best, this simple 
procedure is physically sound; however, it is not 
computationally robust if the body imparts a great 
amount of disturbance to the free stream, such as 
in a hypersonic stream. Hence, the code uses an 
alternate procedure that prescribes the essential 
features of the inviscid flowfield and starts the 
calculation with a rather crude initial approxi- 
mation. 

To accommodate the various shock shapes and 
body configurations displayed in Fig. 1, a reason- 
ably accurate procedure for inviscid flow calcula- 
tions has been implemented into the code. The 
stagnation properties and jocation are determined 
on the pitch plane from n-V = 0, then, taking 
into account the angle of attack, a body angle is 
found to be equal to 


0 = cos cos (a - a^) + n-k sin (a - a^) 


With 0, a distribution of Mach number is estimated 
to be M = |2 sin 0| or M = 1 + |sin 0( for 0 a 
n/4. Using the Mach distribution and the stag- 
nation properties, isentropic expansion formulas 
can be applied to calculate the remaining flow 
variable on the body. This procedure is found to 
be more versatile than the one based on the Newto- 
nian pressure and the conservation law of total 
enthalpy on the wall, e.g. 


+ ‘>•6'^' = *„ + 0.6(V, - n(n-vjf 
P„ = YP„/(y - 


The flow vector tangent to the wall is readjusted 
in accordance with Eq. (6) discussed in sec. 2.4. 
The shock standoff distance at the stagnation 
point is obtained simply from % = 0A5M^/{M^ - 
1)11; 5 ^ is multiplied by 0.5 if an equilibrium- 
air model is used. 

The complete shape is estimated from S = B + 
8o exp(a|^ + a cos <|) for Fig. 1a, 5 = B + 8o( 1 + 
a(2/n(^ - ar cos 0)2) for Fig Ic, and S = B + 8o(l 
+ a cos 0 for Fig. 1b, where a is an input param- 
eter. The shock location takes the following form 
for Af ^ 1 , S = d/(cos ^ + 1/c), a quadratic curve 
intersecting ^ = 0 and n/2 lines; d and c are pa- 
rameters to be determined for each case consid- 
ered. Once the shock is known, the stationary RH 
equations are used to calculate shock variables 
which are used with the wall variables to deter- 
mine the variables inside the shock layer by 
linear interpolations. 

For viscous flow calculations, it is recom- 
mended to start from a nearly converged inviscid 
calculation for which the major portion of the 
shock layer is well established. In this way, a 
substantial reduction of computation time can be 
accomplished. 
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2.6 Wall Shear. Heat Flux, and L/D Calculations 

Vector components of velocity, stress, and 
heat flux can be easily obtained in the wall in- 
trinsic frame from those in the spherical frame. 
With reference to Table 1, a transformation matrix 
T may be defined to relate the components in both 
frames by 



Likewise, a tensor can be transformed accordingly 
by 


n' = rnr^ (8) 


where the prime denotes the intrinsic frame and 
superscript T denotes transpose of the transforma- 
tion matrix. The expanded expression of stream- 
wise and crossflowwise shear stress components and 
of the body-normal heat-flux component are as 
follows: 


If a and ar are nonzero, the drag and lift forces 
can be obtained in the wind frame from 


L = L cos(a - sin (a - a^) 

D = sin (a — + D cos (a — a^) 

Implicit Difference Method 

3.1 Difference Approximations 

The computation space (jc,y, 2 > is discretized 
into a three-dimensional network of grid points, 
which are defined by 


= (m - 2)Ax, y^ = in- 2)Ay, z^ = (e~ 2)Az ( 10) 


with Ax = Xtnax/(f^c - 2), Ay = 1/(nc - 2), Az = 
n/(/c - 2), and m= 1, 2, ... me, « = 2, 3» ... 
nc, and f = 1, ... /cp, where /cp = /c + 1. The 
subscripts m and € cover a range of grid points 
beyond the computation domain (Fig. 2) in order to 
incorporate the boundary conditions. Equation 
(10) indicates that any point in the computation 
space can be addressed by the intersection of 
three independent families of coordinate lines, 
but the metrics and the dependent variables still 
require three-dimensional arrays such as pn,m,€ = 

The partial derivatives are approximated by 
central- and one-sided-difference formulas 


V = (Pm + l-Pm-l)^2Ax 


n , = + Ln. ) + nJt.n . 

nt 2 2' 1 




Cp = (p™.i 


p )/Ax, 8 p = (p - p -)/Ax 


"ni = + ^2% + W + '’2"nn 

+ 

% = + Vn + 


The drag and lift are obtained by integrating the 
pressure and shear stress on the body surface from 


D = \ f "^’‘{n-Kn + t-Kn , + b-Kn | 

Jo Jo \ "" 


(9) 


- = I (n-Ia 


+ Mn + b ln . 
nn nt no 




where da = sin ^ dl, and the scalar products 
are, for example, 


wherein subscripts n and € are excluded for brev- 
ity. The one-sided quotients and subscripts m ± 
1/2 are used primarily for the diffusive flux 
terms in Eqs. (1) and (2). For example, 


where 






ti'K = sin ^ — ^2 oos ^ 

7vl — njCos^oos^+ n^sin^cos^— n^^sin^ 


Special expressions for diffusive terms are needed 
on X = 0 wherein ^ 3 = 0 , but details will not be 
presented here. 
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^,2 ADI Factorization Technique 

The resultant finite-difference counterpart 
of Eq. (1) is 


6/u + 6 /+6g + 6 /» + r , = 0 
t X z n,m,f 


( 11 ) 


where 


8,-*- « = (u* - u* ,)!M + 0(A<2) 

t n,m,t 


8 f= , ,)/2 A*, etc. 


The spatial approximations are solved simul- 

iit+i 

taneousLy with the unknown vector Un,m,e to allow 
greater changes of and, hence, to result in a 
faster rate of convergence from initial approxima- 
tions to the final solution. A factorization 
technique due to Beam and Warming (8) which empha- 
sizes the noniterative nature of the well- 
established ADI method and solves for the incre- 
mental vector of the conservative variable u has 
been adopted and modified. In this version, 
because the primitive variables v - (p,a,u,«;,c)^ 
are used, the Jacobian matrices have simpler 
expression and involve less arithmetic. The 
solution procedure comprises five steps. 


It Is seen in Eq. (12) that the first step is 
to find the residual explicitly, then to smooth or 

k 

to filter out the localized errors in by 

subsequent implicit steps performed along each 
coordinate. Obviously, the implicit calculations 
are not needed for numerical stability considera- 
tion alone if meet the Couran t -Fried rich- 

Lewy (CFL) criterion, or any step can be deleted 
if > CN(Mn,m,€)minf Where i = ft, m, or f 

and CN denotes the Courant number. Nevertheless, 
the primary function of implicit solutions seems 
to be more in the acceleration of the overall rate 
of convergence than in maintaining solution 
stability. Numerical experiments have suggested 
that implicit solution is even preferable (without 
concern in computation cost) to the explicit 
solution because the effect of boundaries can be 
transmitted across the entire line. To illustrate 
the manner in which boundary conditions are incor- 
porated into the tridiagonal system of equations, 
the third equation of (12) is rewritten in alge- 
braic form for index n. 




a Av** , + 6 Au + c Ad . = * 

n n-l n n n n+1 


a Ad . + 6 Ad = 4> 

ncm ncm—l ncm ncm nctn 


(13) 


where 


[rl + 8 A - 8JD + =xP-‘(A«* „,) 

(t/ + 6B - 8JE + = xAr^^^ 

(t1 + 8C - 8JF + trf,7))Ar* + ‘^ = xAr”^^ 


n,m,c «,m,c 


( 12 ) 


a = — fs , + — (B , +i«(,)) 
^ 2 Ay V ” “ ^ Ay ^ ~ ^ ^ ' 


6 = X + -ziE + xd.) 

n n / A \2 « ^ 

(Ay) 


= —(b — (B 4., + 

n 2Ay V Av ' 


Ay 


b =6 for tt.D, ID 

ncm ncm 


= b +c T IT forp 

ncm ncm ncm nc 


where x is the inverse of the local Af multiplied 
by a constant Courant number, and dj and are the 
implicit and explicit damping coefficients to be 
discussed later. Jacobian matrices P, A, B, C, 

O, Ef and F have analytical expressions. Each of 
the three equations in the middle of Eq. (12) 
represents a tridiagonal system of linear equa- 
tions with 5 by 5 matrix coefficients. A standard 
algorithm has been used earlier in Ref. 9. 


Equation (13) solves for Ad„ from n • 2 
through ncm = fic-1; on the wall, Au„c = 

= Aencm = 0 no-slip and isothermal conditions 
and Pnc = Pncm implies Ap„c. It is noteworthy to 
point out that implementation of boundary condi- 
tions is relatively inconvenient for conservative 
variables. 
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3.3 Jacobian Matrices 

The nonlinear difference equation (11) is 
linearized according to Newton’s method prior to 
the factoring of spatial operators into Eq. (12). 
The iteration steps of solving a set of nonlinear 
equations can be minimized if the Jacobian matri- 
ces are exact and updated at each step. By def- 


inition, Au = P Ay, A/ = PA Ay, Ag = PB Ay, A/r = 
PC Ay, A/*y = PD A(8y*), Ag„ = PE A(5y^), Ah^ = PF 
A{^Vz)j where subscript y refers to the diffusive 
fluxes, and x, y, and z denote the respective 
derivatives. Note that the linearized difference 
equation contains no Jacobian matrices due to 
mixed derivatives, because of the nature of the 
ADI factorization technique. 
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where ^ = y * \ z ] ^ p> {ge) for both ideal-gas and equilibrium-air models. 


8 



3.4 Numerical Damping and Time Increment 

Two types of damping terms are included in 
Eq. (12); a second-order term with constant value 
is used in the implicit operators, whereas a mix- 
ture of both second- and fourth-order terms is 
used in the explicit operators. The implicit 
damping terms are of the following form, 

- d ,(8 +8 +8 

/ XX yy zz 


the diagonal dominance of the matrix. The allow- 
able ‘CN is also determined in part by the close- 
ness of initial approximations to the solution; a 
smaller value is often needed at the beginning of 
the iterations, then a larger value is used later 
to optimize the convergence rate. If the maximum 
incremental variable of Ap/p is used to gauge CN 
for the next iteration, then is heuristically 
controlled by the outcome of each iteration. 
Therefore 


Thus, this value enhances the diagonal dominance 
of the matrix and effectively extends the range of 
CN numbers. Although there is discernible lag of 
convergence history initially when the implicit 
damping is used, the convergence rate is faster as 
the solution reaches the asymptotic state. The 
best aspect of all is that the final results are 
nearly independent of dj. 

The use of explicit damping terms, however, 
exerts a profound effect on the solution stability 
and accuracy, because the added numerical dissipa- 
tions become a permanent part of Eq. (12). The 
positive aspect is to smooth out spurious oscilla- 
tion of results caused by the indiscriminating 
difference across flow gradients; the negative 
shortcomings are that they tend to mingle with the 
physical dissipation and markedly change the nat- 
ural phenomena. For the blunt-body flowfield 
computation, the explicit damping terms can be 
selectively used and are beneficial rather than 
pernicious to the final solution. The second- 
order terms are only used for the inviscid portion 
of the shock layer. They consist of an expression 
giving a blanket smoothing and another expression 
depending on the pressure gradient, such as 


Ai = CNm . /(I + e) 


where e = max (Ap/p)„,„i^^. The local-time increment 
for a constant CN 


U 


CN (Aq, q A^ q sin ^ AO 


n,m,€ 




1 / 2 ] 


(15) 


is to be used for all grid points in Eq. (12). 
However, the shock position is still updated by 


gk+i _ gk 






( 16 ) 


The shock speed is smoothed before entering in Eq. 
( 16 ) by 




+ d (8 


+ 8 + 8 (8 
yy zz p x: 


+ 8 +8 
yy 2 




where, for y < C 3 , c?e ,2 = 0.02 to 0 . 2 , dp = 

\hx‘^pK^p)\ ^ and 03 is a parameter delineating the 
boundary layer from the inviscid shock-layer flow. 
The fourth-order damping terms are used for the 
entire domain excluding the points on and adjacent 
to the boundaries. They are of the following 
form. 


-d ( 6 ‘‘ + 8 ‘‘ + 8 ‘bu* 

ex y z 

where dg ^4 - min (Ai,de^). 

The time increment Ai = 1/x in Eqs. (11) and 
(12) is a multiple of (AOmm given by the CFL 
conditions for stability. 


A^ = CN{M . ) = 


CN min (Aq, q A^, q sin ^ AQ 
min 1.5 (a + . 


(14) 


where a is the speed of sound. The Courant number 
CN is greater than unity but its magnitude 
depends partly on the values of dj and rfg. It is 
seen that the fourth- and second-order damping on 
the right-hand side are independent of but the 
damping on the left-hand side as well as x affects 


if ^ + a cos ^ > Cl, where ci ranges from 0.5 to 
0.9. 


4. Code Capabilities and Limitations 

As indicated in sec. 1, the code is developed 
to provide flowfield data for the body configura- 
tions shown in Fig. 1. The body contour may be 
given either analytically or by a lofting proce- 
dure. The contour is allowed to be indented, if a 
uniform x,n is considered adequate to resolve the 
mild changes of body curvature. A choice can be 
made between a perfect gas and equilibrium air, 
but the turbulent eddy viscosity has not been 
incorporated. The computation domain is prefer- 
ably no greater than 45"" in the angular span of 
the cone generator since the shock fitting becomes 
less accurate as the angle cos”i (i-Voo/|Va>|) de- 
creases. The speed of the incoming flow is 
unrestricted for the perfect-gas solution; how- 
ever, it may impede the accuracy of shock fitting 
because, with an additional shock speed intro- 
duced, the resultant temperature behind the shock 
could be higher than the range of the curve fit 
for equilibrium-air properties. The location of 
the coordinate origin also controls the size of 
domain which is required to cover the entire 
subsonic region inside the shock layer. For 
instance, with an angle of attack of 40®, the 
sonic line is found to extend as far as ZRn from 
the nosetip. Higher flow incidence angle may 
cause the grid lines to intersect the shock and 
the body at excessively skewed angles and to give 
rise to computational errors. The configuration 
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of Fig. la represents the nose of a high lift-to- 
drag (L/D) vehicle; hence, both nose and axial- 
plane data are needed as shown in Ref. 9. The 
configuration of Fig. lb, suitable to start a 
method-of-characteristics solution, was considered 
in Ref, 2. The configuration illustrated by Fig. 
1c has sonic shoulder of varying degrees of 
curvature and a wide-angle conic frontal surface. 
This type of body is being investigated for 
aerobraking near the upper Earth atmosphere. The 
following discussion will highlight the parameters 
used to set up the computation. 

The Euler code is robust and accurate while 
forgiving to the spatial resolution, CN number, 
and damping coefficients. As long as the outflow 
boundary is well outside the subsonic zone, coarse 
grids of 5 by 12 by 0 and 10 by 12 by 7 for 2D and 
3D cases, respectively, are sufficient to yield 
results accurate to the third digit. The calcula- 
tion is stable with CN = 3 to 15 for global-time" 
increment and CN = 3 to 6 for local-time incre- 
ment. The time steps are between k = 100 to 300 
to achieve convergence. The criteria of solution 
convergences are based on the maximum incremental 
variable c = max |Ap/p|, the maximum range of shock 
speed St = |max St - min and the maximum 
deviation of total enthalpy A//. For coarse grid 
solution, 0.001 ^ z ^ 0.01, 0.01 <: ^ 0.1, and 

A// ^ 0.005 are considered sufficient. The 
criteria can be lowered according to the reduced 
grid spacing. The seemingly high threshold for St 
is due to the inaccuracy of shock fitting near the 
upper end of the shock. It should be noted that 
the corresponding root-mean-square (rms) value of 
Si is generally one order of magnitude less than 
5^. (15) For most calculations, the time increment 
and damping coefficients are observed to exert 
noticeable effect on the converged results; it is 
a recognized characteristic associated with the 
ADI technique. The explicit damping terms enhance 
the computational stability as well as the code 
robustness. Both second- and fourth-order damping 
are used and an additional shock-speed-smoothing 
scheme is introduced to ensure that the shock and 
flow variables have no anomalous oscillation. The 
damping coefficient de ^ 0.1 is recommended for 
inviscid calculations. The effect of damping 
helps to expel the errors outside the domain, as 
Aif may vary from 0.2 to 0.005 during the time 
integration. If AH has a constant level of 
error, a slightly greater de may be needed. 

The viscous calculations can be performed 
from the results of their inviscid counterparts 
or, alternatively, by going through a straight- 
forward coarse-to-fine grid sequence. Both have 
been used and found to be cost-effective. After 
interpolating the coarse inviscid grid to grids of 
28 by 10 by 0 or 28 by 12 by 7 on a nonuniform 
grid exponentially stretched toward the shock, the 
grid spacing adjacent to the wall is estimated by 
Ay = 0 . 2 /?at/ ( P i?e ) 1/2 ^ where p is a parameter equal 
to the subscripts st and w refer to the 

stagnation and the wall, respectively. As more 
grid points are packed toward the wall, sparse 
grid points are seen near the shock. To prevent 
shock-fitting errors caused by the stretched grid 
from accumulating, more than 28 points may be 
occasionally needed to maintain St to a negligible 
level. The grid spacing Ax and Aa is not as 
important to the accuracy unless around a sonic 
shoulder as in Fig. 1c or for the indented nose 
investigated in Ref. 10. One should be cautioned 
that the placement of the outflow line is more 


critical than for the inviscid calculation; 
placement is preferably on the axial plane, where 
flow is less accelerating than on other parts of 
the body. The local-time increment is advanta- 
geous to use since St is smaller while the rate of 
convergence for the heat-transfer and friction 
coefficients, Ch and C>, respectively, is faster. 
The Courant number is in the range of CN = 20 to 
35, a factor of 2 lower than that used for global - 
time increment. The criterion for convergence is 
the same as that for inviscid computations plus 
the check made to the values of C// and C>. 

Although the 3D NS code has performed excep- 
tionally well in terms of the grid points and 
integration steps used for the strongly asymmetric 
flowfield considered, surprises do occur in lack 
of convergence, solution smoothness, and accuracy 
of results. The remedies to overcome these diffi- 
culties are 1) reducing grid spacing and CH, 2) 
increasing t/e» and 3) selecting a different 
The following section will demonstrate much of the 
usefulness and the versatility of the code for a 
wide range of practical problems. 

5. Discussion of Sample Cases 

Five problems have been selected to represent 
a partial spectrum of current interests; their 
free-stream conditions, computation domain, grid, 
and body shape are given in Table 2 in accordance 
with the required computational effort. All are 
laminar flow and isothermal wall, and the computa- 
tion domain contains an axial plane at Zo. Coarse 
grids were used so that a scalar computer such as 
Univac 1182 will have sufficient size and speed to 
execute simple 2D cases in time-sharing mode and 
3D cases overnight. Before addressing each indi- 
vidual case, we will explore some basic issues 
that have revealed interesting and noteworthy 
observations, which are shared by the shock- 
fitting implicit numerical procedure for both 
Euler and Navier-Stokes codes. 

The convergence rate, as expected, varies 
with the time increment even though the shock 
position is decoupled from the implicit calcula- 
tion. Shown in Fig. 4a is the convergence history 
of the time integration for an inviscid flow over 
a sphere at a = 0*^ and M = 5.94 obtained with CN 
=3, 9, and 15. This Naval Surface Weapon Center 
(NSWC) case is designated NSWC-1. Using global- 
time increment, the higher value of CN leads to 
lower value of 5^ and e over nearly all time steps 
except near k = 200. At the end of computation, 
both and e are seen to level off where the 
errors due to truncation, linearization, factor- 
ization, and numerical damping have reached equi- 
librium. Further calculations did not alter the 
level appreciably. The convergence history 
obtained from the calculations using local-time 
increment is plotted in Fig. 4b. Despite the 
lower CN = 6 value permitted by local-time incre- 
ment, the rate of convergence is seen as rapidly 
decreasing as that in Fig. 4a for CH = 15 using 
global time. As marked in Figs. 4a and 4b, the 
shock standoff distance and the pressure on the 
stagnation point only differ by from each 
other. Both cases (8 by 12) are considered very 
accurate for AH of less than 0.1?, Discernible 
advantages of accelerating the convergence process 
are noted for viscous flow computations. Figure 
4c shows the convergence rate for both local- and 
global-time integration. The value of St is also 
smaller in local-time calculation. The conver- 



gence rate for global-time integration does not 
speed up as much as expected for CN = 50. Hence, 
local- time integration has been used for all 
results discussed later. The converged Cjj value 
of 0.0042 compared favorably with the value 0.0043 
obtained by Hsieh.(15) 

The significance of using local time as 
opposed to global time is further exemplified by 
Fig. 5, wherein the shock speeds are compared. 

The local -time solution gives a uniformly lower 
shock speed than those obtained with global-time 
calculation. The greatest St occurs at ^ = n/2 
for both cases. The error in S<, estimated by 
|5i|/(\/YAfa) , imposes little upstream influence 
because of its proximity to the outflow line. 

It is inherently associated with the finite- 
difference shock-fitting scheme that is an approx- 
imation to the original method-of-characteristics 
shock fitting advocated by Moretti.(2) Regardless 
of the error magnitude, the shock locations are 
very close to those shown in Ref. 1 for Af = 3: 

8„ = 0.229 vs. 0.23 and 8„/2 = 1.04 vs. 0.9. 

Another important observation is that the con- 
verged results obtained using different damping 
coefficient or CN values may not be in perfect 
agreement. 

A comparison of Ch and Cp distributions be- 
tween the NS and the boundary-layer (BL) results 
is shown in Fig. 6 for the Arnold Engineering 
Development Center (AEDC) case at a = 0®. The 
inviscid and viscous pressure distributions are 
almost identical for this Re range. The NS code 
has predicted higher C> over the entire wall, but 
yielded a C// distribution that coincided with 
that of the BL code(l6) except at the stagnation 
point. The outflow boundary = 98® is obtained 
by means of two-point extrapolation. No adverse 
upstream influence is noticed since the flow no 
longer moves as fast as that on the spherical 
nose. 

The preceding examples were obtained with 
= 0 and deA ^ 0. If d^^z > 0.0, St and c would 
have shown lower magnitude, but the calculation 
was stable without using For three- 

dimensional cases, however, de;z is important in 
stabilizing the calculation. To examine its 
positive effects, the AEDC case was computed using 
three different values and the history of 
convergence is presented in Fig. 7. An inviscid 
calculation was obtained with dgz = 0.1 from k = 0 
to 200, then the viscous calculation proceeded at 
different dg ,2 values. The influence of dg ^2 to the 
stability is self-explanatory. Nevertheless, the 
explicit damping is Incapable of removing the 
spurious oscillations along a coordinate line. 

For instance, the C// distributions were noted to 
swing up and down during the calculation with CN 
= 50 but to be smooth with CN = 30. The results 
at k - 600 (including 200 steps for inviscid 
calculation) are plotted in Fig. 8, but not shown 
is that the shock standoff distance near 0 = 10® 
also oscillates with CN = 50. The value of Ch 
depends on the temperature gradient at the wall 
and thus is directly proportional to the shock- 
layer thickness at that point. Therefore, the 
accuracy of the viscous parameter, being very 
sensitive to the inviscid layer and shock calcu- 
lations, should be used to evaluate the quality of 
the overall solutions. In contrast, the wall 
pressure distribution is relatively insensitive to 
shock shape and diffusive fluxes and cannot be 
used for that purpose. 


The NSWC case 2 (Figs. 9 and 10) was consid- 
ered to certify the conventional approach of 
finding starting data for aft-body calculation at 
Q > 0®. Since the flowfield over a sphere-cone is 
symmetric to the wind axis, an interpolation of 
the available 2D axisymmetric results will provide 
data on an axial plane for 3D cases. By selecting 
Zo = -0.5, the outflow line at ^ = n/2 is quite 
close to the sonic line on the windward plane and, 
thus, may present upstream propagation of extrapo- 
lation error at the downstream boundary. The one- 
point extrapolation is deemed more stable than the 
two-point extrapolation that is used for the 
supersonic region. The pressure values are dif- 
ferent at the stagnation point between the invis- 
cid and viscous calculations. In Fig. 10, the 
heat-transfer distribution exhibits a small 
oscillation caused by the shock locations. More 
oscillations appeared in Ch and St when a grid of 
28 by 12 by 3 was used instead of 33 by 12 by 3. 
Five additional grid points between the shock 
layer had resulted in more accurate shock fitting. 
Favorable comparison of the axial-plane data is 
found for 2D and 3D calculations at a = 10®. 

The AEDC case (Figs. 11 and 12) refers to the 
test environment of tunnel B, which was used 
extensively to develop an aerodynamic data base 
for the Space Shuttle Orbiter. For a simple 
hemispherical cone, the axial plane was placed at 
Zo = 0.5 to minimize the upstream influence within 
the boundary layer. The Reynolds number is not 
high, and the wall temperature is about half the 
total temperature; hence, a grid of 28 by 12 by 5 
is adequate to fit the shock and to resolve the 
crossflow gradient. In Fig. 12, a mild hump is 
displayed on the Ch distribution that can be 
attributed to the shock shape near ^ = 0. A few 
more grid points In q may smooth out both. The 
inviscid and viscous results are otherwise con- 
sidered excellent. The data generated on the 
axial plane at = 0.5 were used to initialize 
both the Euler and the parabolized NS computation 
by the author^ to an axial station equal to Z© = 
30. 

The geometry of the Space Shuttle Orbiter was 
defined by a lofting technique using blueprint 
configuration. The Orbiter nose consists of a 
spherical cap at the noset ip and nonanalytical 
shape. The nonspherical configuration and a s 40® 
necessitate 3D NS calculations not only to deter- 
mine Ch values but to find initial data for 
parabolized NS calculations. Under equilibrium- 
air assumption, the subsonic zone is well within 
the domain located at Z^ = 0.5. Using a grid of 8 
by 16 by 8 and after 300 steps, the inviscid 
results have St = 2.9 and e = 0.009. The viscous 
calculations continued with a grid of 28 by 16 by 
8 and another 300 steps of integration. The 
results have St = 0.4 and Ch = 0.04. The Mach 
contours of the shock-layer flow are shown in Fig. 
13. There are marked differences between inviscid 
and viscous results which may be attributed to the 
vortical flow enhanced by the presence of a vis- 
cous layer on the leeward side of the body. 

Another interesting pattern is that the locations 
of maximum p and maximum Ch are not necessarily at 
the same point on the body. (See Fig. 13.) 

The grid points used to approximate the body 
and the shock are illustrated by Fig. 14 in two 
perspective views. They are very coarse for the 
size of the shock layer. The normalized data for 
the fifth Space Transportation System flight (STS- 
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5) are also presented. The discrepancy between 
the NS prediction and the data was also noted 
between the BL results and the data. No pressure 
data are available for comparison; however, the 
viscous pressure is somewhat higher than inviscid 
pressure on the windward side. 

The aerobraking orbital transfer vehicle 
(AOTV) configuration represents the front portion 
of an ellipsoid/60'' cone proposed by the NASA 
Lyndon B. Johnson Space Center for a flight 
experiment. The complete body consists of the 
front portion, a toroidal shoulder on a 15° raked 
plane, and a payload or propulsion bay. In the 
flowfield computation, the toroidal shoulder is 
replaced by a circular shoulder to facilitate 
computation. The shock and body grid network is 
shown in Fig. 15, together with the shock layer on 
the pitch plane. The purpose of Fig, 16 is to 
point out that the heat flux to the front surface 
is strongly affected by the discontinuities in 
body curvature, whereas the pressure is much less 
sensitive. The substantial differences between 
the inviscid and viscous pressure distributions 
(Fig. 16) imply that at such high altitudes, the 
viscous effects to the aerodynamic characteristics 
of the vehicle cannot be ignored. The computation 
time is 256 minutes, which accounts for both 200 
and 300 iterations, respectively, for inviscid and 
viscous grids. 

Concluding Remarks 

Since the reports published in 1974 and 1975 
of the Navier-Stokes solutions for asymmetric 
entry flow over blunt bodies, significant improve- 
ments have been made to the original 3D NS code. 
Although the formulation is still cast in curvi- 
linear-orthogonal coordinates and uses a special 
equation to solve for the singular axis, the ex- 
plicit scheme has been replaced by an ADI factor- 
ization technique combined with numerical damping 
and local-time increment. The code has a reliable 
and cost-saving procedure stepping from coarse to 
fine grids for solving viscous flow. Not only has 
the computational effort been reduced by a factor 
of 100 or more, the solution accuracy is also 
improved because of better understanding of the 
methodology. This code has a formulation, a 
numerical method, and an algebraic grid that are 
consistent to the level of accuracy desired for 
evaluating engineering design of space vehicles 
experiencing peak heating at high altitudes. We 
have finally succeeded in the stratagem by first 
developing an accurate formulation and then con- 
tinuously improving the computational efficiency 
as better methods become available. 
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Appendix A ~ Derivation of Ea. (2) 


Equation (1) is undeterminate or singular 
along the Z-axis because of sin ^ = 0 and w = 0. 
^n the removal of sin ^ or ^ from the flux terms, 
has a complicated expression and needs special 
analysis. Consider a typical component of S; the 
mixed derivative of second order results in 


(pw^)^ = ( A1) 

But the velocity components are related on the Z- 
axis by 


since 

{ )^ = ( );c + )y» ( )< = ( )z Jcf, )y» 

( )n = ( ). 

The final form is 


u = u(0,n,0) cos^,v- v(0,n,0),«; = -tt(0,n,0) sin < 

Hence, when evaluating Eq. (A1) on the pitch plane 
< = 0, one observes that 

«; = 0, = -u(0»q»0), - 0, = -a^(0,q,0) 


This term, in conjunction with (f + r)/sin ^ from 
Eq. (V), completes the derivation of the convec- 
tive portion of Eq. (2). 

An evaluation of the diffusive fluxes on the 
Z-axis is not pursued analytically; instead, an 
approximated expression is derived using differ- 
ence formula and grid notation presented in sec. 
3.1. Let H be the diffusive fluxes. Then 


and “ -p^w(0,n.0) - pu^(0,n,0) -(pu*)^ 


zx zy 


where u* - u(0,q,0). Thus, the convective h flux 
becomes 


where 




ti(o,n,o) - 


(0,11,0) u 


” pu “ 


“ P 

pu2 

, u - 

pu 

puu 


pu 

0 


0 

^ pwe + 


.Pe +p_ 


» , = H ' ("„ A 1 + ]/(2 ^ 

“■ (^n+1,24 ^ 
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Table 1 Unit vectors in reference coordinate frames 



Wind frame 

Cartesian frame 

Spherical frame 


U.K 

I.J.K 

ij,k 

Cartesian 

Icosa + Ksina 

1 

i cos ^ cos ^ + j sin ^ cos ^ -- A sin ^ 

frame 





J 

1 

i cos ^ sin ^ + j sin ^ sin ^ + 6 cos ^ 


-/ sin a + K cos a 

1 

isin^-jcos^ 



7 cos ^ cos ^ + J cos ^ sin < + 7C sin ^ 

■ ■ ^ 

Spherical 




frame 


7 sin ^ cos ^ + J sin t sin C, -K cos ^ 

1 

ij,k^ 






-7 sin ^ + J cos ^ 

1 


ini-hjn^-\- Ana 

Intrinsic 

frame 

’ ’ 161+^62 + ^63 


e.g., 

V=iu+ jv + kw = lU + JV + KWr 

U^IV=^Iiu + I‘jv + Ikw 
Table 2 List of sample cases 




Case 

M 

a, deg 

Qr, deg 

Zd“ 

Grid 

Modelb 

Tw 


ft 

NSWC-l 

2.97 

0 

0 

0 

8byl2 

P 

1 

1.01 X106 

0.1823 


5.94 




28byl2 





NSWC-2 

5.94 

10 

0 

-0.6 

8byl2by3 

P 

4.4 

1.01 X 106 

0.1823 






33byl2by3 





AEDC 

8 

30 

0 

0.5 

8byl2by5 

P 

6.5 

0.71 X 106 

0.0413 






28byl2by5 





Orbiter 

22 

40.8 

0 

0.5 

8byl6by8 

E 

8.5 

14 000 

2.36 






28 by 16 by 8 





AOTV 

34.8 

0 

15 

NA 

10 by 18 by 5 

P 

8.4 

8272 

10 






28 by 18 by 5 






aLocation of the Cartesian frame relative to the center of the nose curvature on Z-axis. 
bP-perfect gas (y = 1.4), E - equilibrium air. 
cFree-stream Reynolds number per foot. 

4Nose radius on Z-axis. 
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NSWC-1: Ms 5.94, 

a sO®, d2 = 0> d 4 «p = 0.1 

8 K 12 K 0 

CN 

^0 Po 

3 

0.1489 45.81 

L 9 

0.1474 45.54 

\ 16 

0.1473 45.95 

\ GLOBALt 


NSWC-1 

M s 6.94. a = 0® 

6 » 12 » 0 LOCAL t CN » 6 
d = 0.1 
6q = 0,1478 
K/2 = 0.647 


NSWC-1 
M - 6.94. a - 0® 


Oq ^7t/2 
. LOCAL t 0.1502 0.609 
• GLOBAL t 0.1519 0,613 


Fig. 4 Convergence history for 2D NSWC-1 case. 

a) Effects of time -increment step size 
on inviscid solution, b) Effects of 
local-time integration on inviscid 
solution, c) Comparison between 
global- and local-time convergence 
rates for viscous case. 





PRESSURE 




-1.0 -0.5 0 .5 1.0 -1.0 -0.5 0 .5 1.0 

(2t)max 


. 5 Inviscid flowfield results for the NSWC-1 
case (M = 2.97, a = 0°); = 

0.58102. Shock speed: 0.003187 (left 

plot); 0.005339 (right plot). 


M = 8, R|m = 0.0413 FT. T,* = 99« R, Re/FT = 707 000 
A8=8», Ar = s/6 

NS EULER O BL 


e, DEG 


b) 

(?. DEG 

c) 

e, DEG 

Comparison of NS and 
tion coefficient. 

BL 

c) 

results for an 
Heat-transfer 

AEDC 2D case, 
coefficient. 

a) Wall pressure 

distribution, b) Fric 






e, DEG 

Fig. 8 Comparison of friction and heat-transfer 
coefficients for an AEDC 3D case using 
different time-step sizes. 


GRID 


MACH 


PRESSURE 








Fig. 9 Inviscid flowfield results for the NSWC-2 
case (M = 5.94, a = 10°); Xmax = 

0.9496. Shock speed, 0.017368 (left 
plot); stagnation pressure, 44.5016 
(right plot). 
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Fig. 10 Viscous flowfield results for the NSWC-2 
case {M - 5.94, a = 10°); X„,ax = 
0.9496. Stagnation Stanton number, 
0.00374 (left plot); stagnation 
pressure, 43.2118 (right plot). 


Fig. 11 Inviscid flowfield results for the AEDC 
case (M = 8, a = 30°) » ~ 1«2161; 

stagnation pressure, 77.2282. 



GRID PRESSURE MACH 



Fig. 12 Viscous flowfield results for the AEDC 
case (Af = 8, a = 30°). X^ax = 1.2161; 
stagnation pressure, 82.1822 (right 
plot); stagnation Stanton number, 
0,03695 (left plot). 



Fig. 13 Inviscid and viscous flowfield results 
for the Orbiter case (M = 22, 
a = 40.8°, equilibrium air); Xmax ~ 
1.613. Stagnation Stanton number, 
0.03796 (left plot); stagnation 
pressure, 663.748 (right plot). 
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Fig. 14 Shock shape and comparison with flight 

data for the Orb iter case (a = 40.8®, t 
= 650 sec). Xmax = 1.613; stagnation 
Stanton number, 0.040736. 



Fig. 15 Shock shape for the OTV case (M = 34.8, 
07.= 15®). 



P 



PRESSURE 

CONTOURS 




Fig. 16 Inviscid and viscous flowfield results 
for the OTV case; X„^ax = 2.493. 
Stagnation Stanton number, 0.02512 
(lower right); stagnation pressure: 
1540.56 Cupper left); 1553.33 (lower 
left, viscous case). 
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